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PACS 05 . 45 . Tp - Time series analysis 

PACS 05 . 45 . Xt - Synchronization; coupled oscillators 

PACS 89 . 75 . He - Networks and genealogical trees 

Abstract. - Synchronisation between coupled oscillatory systems is a common phenomenon 
in many natural as well as technical systems. Varying the strength of coupling often leads to 
qualitative changes in the complex dynamics of the mutually coupled systems including different 
types of synchronisation such as phase, lag, generalised, or even complete synchronisation. Here, 
we study the geometric signatures of coupling along with the onset of generalised synchronisation 
between two coupled chaotic oscillators by mapping the systems' individual as well as joint 
recurrences in phase space to a complex network. For a paradigmatic continuous-time model 
system, the transitivity properties of the resulting joint recurrence networks display distinct 
variations associated with changes in the structural similarity between different parts of the 
considered trajectories. They therefore provide a useful indicator for the emergence of generalised 
synchronisation. 

This paper is decidated to the 25th anniversary of the introduction of recurrence plots by 
Eckmann et al. (EPL, 4 (1987), 973). 
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Introduction. — The scientific interest in studying 
and qualitatively understanding synchronisation between 
coupled oscillatory systems traces back at least to observa- 
tions of weakly coupled pendulum clocks made by Christi- 
aan Huygens in the ITth century jl]|2]. In the last decades, 
a detailed quantitative knowledge of synchronisation phe- 
nomena has been established for both periodic and chaotic 
oscillators (3||7]. 

In the study of synchronisation processes between cou- 
pled chaotic oscillators, different phenomena have to be 
distinguished: First, complete synchronisation (CS) [6] 
can occur for identical systems with a sufficiently strong 
coupling. In this case, the trajectories x{t) and y{t) of the 
two coupled systems become identical, i.e. y{t) — x{t). 
Second, generalised synchronisation (GS) [8}{TT] refers to 
a general fixed and deterministic functional relationship 
between both trajectories, y{t) — f{x{t)), where / is a dif- 
feomorphism. This type occurs for non-identical chaotic 
oscillators at rather high coupling strengths. A closely 
related form is lag synchronisation (LS) 12 , where both 



coupled systems evolve in an identical way with a fixed 
mutual time shift as y{t) — x(t — r). Finally, phase syn- 
chronisation (PS) is characterised by a locking between the 
phases of two systems, m(j)y{t) = n(l)x{t) with TO,ri e N, 
whereas the amplitudes may remain uncorrelated. PS is 
typical for non-identical chaotic systems at moderate or 
even weak coupling strengths. 

Even though PS is typically considered a weaker form 
of synchronisation than GS and CS, the sequence of syn- 
chronisation phenomena arising for two coupled oscilla- 
tors with increasing coupling strength depends crucially on 
the structural differences between both systems as well as 
the specific coupling configuration. For several examples, 
it has been demonstrated that as the coupling becomes 
stronger, the dynamics of both systems becomes succes- 
sively synchronised on more and more characteristic time- 
scales 
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This observation implies the emergence of 
PS at relatively low coupling strengths if the characteris- 
tic oscillation frequencies of both systems are sufficiently 
similar, whereas LS, GS, and CS typically occur at higher 
coupling strengths. However, under certain conditions the 
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observed synchronisation scenario can be much more com- 
plex. 

A formal description of GS has been given first only 
for driver-response systems. However, the typical signa- 
tures of GS are not unique to systems with unidirectional 
coupling, but can be observed in a similar way also for 
symmetrically coupled oscillators [l5] (in fact, GS has 
been first described for two systems with bidirectional cou- 
pling [5]). 

For driver-response relationships, the driven system un- 
dergoes systematic changes in its dynamical evolution as 
the coupling strength is increased. These changes can 
be detected by a variety of methods such as the auxil- 
iary system method [9], numerical estimates of the con- 
ditional Lyapunov exponent of the response system as a 
measure for its asymptotic stability To], or by quanti- 
fying the predictability of the driven system. The lat- 
ter class of approaches includes several methods based on 
neighbourhood relationships in phase space (i.e., relying 
on geometric information), with the mutual false near- 
est neighbour method |8 and the synchronisation likeli- 
hood J6| as most prominent examples. However, there 
are specific cases (such as systems exhibiting multistabil- 
ity) where at least some of these methods are not applica- 
ble. Even more, they are not always directly transferable 
to the case of bidirectional coupling configurations. Con- 
sequently, recently several new approaches have been pro- 
posed for studying GS for both uni- and bidirectional cou- 
plings, including methods based on recurrence plots [17| 
or a generalised angle between the trajectories of both sys- 
tems [18]. 

Since synchronisation has distinct effects on the spa- 
tial organisation of driven systems in phase space, geo- 
metric methods are promising and widely applicable can- 
didates for detecting GS from time series data. In this 
work we propose a new geometric method based on a com- 
plex network approach. Specifically, our method is based 
on a graph-theoretical interpretation of joint recurrence 
plots describing the simultaneous occurrence of mutually 
close pairs of state vectors in two or more coupled sys- 
After presenting the details of our approach. 
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tems 

the paradigmatic example of two coupled Rossler oscilla- 
tors in different dynamical regimes are discussed in order 
to demonstrate the performance of the proposed method. 
The obtained results indicate that the complex network- 
based approach has great potentials and may be applied 
to real-world problems in future work. 

Methodology. — In the last decades, there has been 
a rising interest in studying structural properties of net- 
worked systems across disciplines 20 21 . For this pur- 



pose, a variety of measures have been introduced that 
quantitatively characterise different aspects of the ob- 



served networks' complex connectivity properties 22 23 



Besides investigating real-world systems with a physical 
network substrate, recently there have been various suc- 
cessful attempts to applying complex network-based ap- 



proaches also to more general problems of data analysis, 
such as for infering so-called functional networks from mul- 
tivariate time series 24 25 or characterising structural 



properties of time series. For the latter purpose, a mul- 
tiplicity of complementary approaches (such as cycle net- 
works, visibility graphs, and recurrence networks, to men- 
tion only a few) have been proposed [26^30^ , which trans- 
form certain types of time series to the complex network 



domain (see 31 for a corresponding review). 

Among the currently existing types of time series net- 
works, approaches characterising the mutual proximity of 
state vectors in phase space have recently become partic- 
ularly popular, since the resulting networks' structures re- 
flect nontrivial geometric properties of the supposed (low- 
dimensional) attractor underlying the observed dynam- 
ics. A particularly interesting representative of this class 
are recurrence networks [30 32 , which encode the mutual 
proximity of state vectors in phase space arising from dy- 
namical recurrences in the sense of Poincare. Specifically, 
defining neighbourhoods of individual states by consider- 
ing a spatial threshold distance e around each observed 
state vector offers a flexible way for capturing the geo- 
metric backbone of the system under study. Formally, the 
system's finite-time recurrence properties infered from a 
time series {xi}^i are then represented by the binary re- 



currence matrix 



R,,{e) = e(£- 



(1) 



where || • || is the maximum norm in phase space (how- 
ever, other norms could be used here as well), and 0(-) 
denotes the Heaviside function. The visualisation of this 
symmetric matrix is commonly refered to as the recurrence 



plot 33 and has become an important tool for nonlinear 



time series analysis in the last decades |34j. 

Recurrence network (RN) analysis reinterprets matrix 
([I]) as the adjacency matrix of an undirected simple 



graph 30 by setting 



Aij (e) = i?y (e) - 5i. 



(2) 



(where 5ij is Kronecker's delta). The resulting networks' 
properties can be described analytically by interpreting 
RNs as random spatial graphs with the system's invari- 
ant density uniquely determining RN connectivity 35p6 



This solid theoretical foundation allows using RNs as a 
widely applicable tool for detecting dynamical changes in 
non-stationary time series from mathematical models as 
well as real- world applications 37 -39 . One basic, yet im- 



portant characteristic of RNs is the edge density 



P{e) 



1 



N{N -I) 



(3) 



which is a monotonous function of e [40 . In general, a 
suitable resolution of small-scale network features requires 
the choice of a low edge density, typically p < 0.05 36p0 



which still meets the condition p > p{ec) where £c is the 
percolation threshold of the RN [36] . 
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It has been shown that among others, the transitivity 
properties of RNs provide some particularly useful char- 
acteristics measuring the effective local as well as global 



dimensionality of the underlying dynamical system 35 



Since in a given integer-dimensional manifold, the transi- 



tivity properties can be computed analytically 41 , gener- 



alising the corresponding relationship to non-integer val- 
ues allows defining the so-called transitivity dimension 



Dr{e) = 



logr(£) 
log(3/4) 



with the network transitivity |21||22 



Tie) 



(4) 



(5) 



Other global RN characteristics (e.g., average path length 
£ or global clustering coefficient C) have also proven their 
capabilities 30pT , but shall not be further discussed here 
for brevity. Specifically, the dimensionality interpretation 
of T will be the foundation of the approach to studying ge- 
ometric signatures of synchronisation that will be detailed 
below. 

Motivated by the great potentials of RN analysis of indi- 
vidual dynamical systems, a thorough extension to study- 
ing the collective dynamics of two or more coupled systems 
has been recently proposed in terms of inter-system re- 
currence networks (IRNs) [42j. IRNs link the RNs of two 
different systems using the cross-recurrence matrix [43| 



CR,,{e) = e(£ 



(6) 



between the associated time series {xi\ and {uj}. Asym- 
metries in the resulting cross-transitivity measures for in- 
terdependent networks f44' have unveiled geometric sig- 
natures of unidirectional couplings between structurally 
similar systems 42 . Note that IRNs do not require simul- 



taneous observations of the considered systems, but imply 
that they share the same phase space - a condition which 
is often not met in real-world applications. Moreover, the 
application of IRNs is restricted to studying weakly uni- 
directionally coupled systems before the onset of GS 42 



In order to study the geometric signatures accompany- 
ing the onset of GS in some more detail, we suggest util- 
ising another multivariate generalisation of the recurrence 
plot concept, joint recurrence plots 19 , capturing the si- 



multaneous occurrence of recurrences in two dynamical 
systems (Fig. [T]) . The underlying joint recurrence matrix 



Q{£x - \\xi - Xj\\)Q{ey - \\y, - y^W) (7) 



is defined as the point-wise product of the individual sys- 
tems' recurrence matrices. This basic concept can be di- 
rectly generalised to the study oi K > 2 coupled systems. 

Next, as for the traditional RN we reinterpret the joint 
recurrence matrix as the adjacency matrix of a so-called 
joint recurrence network (JRN) by setting 





Figure 1: Schematic illustration of a joint recurrence of two 
dynamical systems X and Y in their respective phase spaces 
PSx and PSy- 



Unlike IRNs, JRNs can be constructed for systems with 
different phase spaces, but explicitly require simultaneous 
observations of all systems and, thus, time series of the 
same length. The system-specific thresholds £k should be 
chosen such that the edge densities pk{ek) of the RNs of 
the individual systems are fixed at the same value p. By 
varying p, a desired edge density pj can be obtained for 
the resulting JRN. 

The properties of JRNs can be widely interpreted in 
a similar way as those of a classical RN in the higher- 
dimensional phase space of the composed system, with the 
exception that spatial proximity is evaluated separately in 
the subspaces belonging to the individual systems. Specifi- 
cally, the transitivity T(£2:, ^y) (cf- Eq. ([5|)ofthe JRN (or, 
for short, joint transitivity Tj) provides a measure for the 
"joint dimensionality" of the composed dynamical system 
(in a general sense). Due to the larger effective dimen- 
sionality of the joint phase space, for two systems X and 
Y in the absence of synchronisation we expect 73s: i Ty S> 
Tj (specifically, = logrj/log(3/4) « Dr^ + Dr^, 

cf. Eq. Q). In turn, when both systems become syn- 
chronised (in the sense of GS), the effective degrees of 
freedom of both systems become mutually locked, result- 
ing in 7j — >■ Tx,Ty, i-e., both systems acting collectively 
as one with the same effective dimensionality as the indi- 
vidual systems. Consequently, studying the ratio between 
the joint transitivity and some mean value (e.g., arith- 
metic or geometric mean) of the individual systems' RN 
transitivities can provide a useful geometric indicator for 
the emergence of GS. Here, we use the arithmetic mean, 
since it is never smaller than the geometric one, implying 
a higher discriminatory power of the transitivity ratio 



Tj 



{Tx+Ty)/2 



(9) 



-^-ij (^2:: ^y) — J Riji^^x i ^y) '^r 



(8) which we will use as an indicator for GS in the following. 
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Results. — As a well-studied paradigmatic model ex- 
hibiting different types of synchronisation phenomena, 
we study the performance of JRNs for two non-identical 
Rossler oscillators [451 



Xi = -(1 + V)X2 - X3 

X2 = (1 + + ax2 + HYxiy2 - X2) 
X3 = b + x^^xi — c) 

m = -(1 - '^)y2 - ys 

2/2 = (1 - i^)yi + ay2 + y'XY{x2 - V2) 

ys = & + 2/3(2/1 -c) 



(10) 



that are diffusively coupled via their second component 
as in j'17 . The detuning parameter v accounts for a fre- 
quency mismatch of the nonlinear oscillations exhibited 
by both systems, making them non-identical and, thus, 
prone to GS. In the following, we will keep it fixed at 
v = 0.02. In turn, for the characteristic parameters a, b 
and c, we will mainly keep the same parameters for both 
systems throughout all further considerations. Specifi- 
cally, we study the geometric signatures of the complex 
synchronisation scenarios arising for two different dynam- 
ical regimes with chaotic oscillations: the phase-coherent 
regime (a = 0.16, & = 0.1, c — 8.5) and the non-phase- 
coherent funnel regime (a = 0.2925,6 — 0.1, c — 8.5). 
In all simulations, we integrate the coupled system with 
step size h = 0.001 until T = 6,000, i.e., for 6,000,000 
time steps. The first 5,000,000 time steps are discarded 
to avoid even very long transients, and the remaining 
data are downsampled by a factor of 200 to obtain RNs 
and JRNs with N — 5,000 vertices. For comparative 
purposes, we also consider the recurrence plot-based in- 
dices CPR and J PR [l7). CPE, an indicator for PS, 
is based on the cross-correlation between the two sys- 
tems' recurrence-based generalised auto-correlation func- 
tions, whereas JPR attempts to identify synchronisation 
by comparing the systems' RN edge densities px.Y to pj, 
which are conjectured to become very similar in case of 
GS. 

Let us first examine a unidirectional coupling configura- 
tion X — > y {pyx ~ 0), where the slightly faster system 
X drives the slower one Y. For two phase-coherent Rossler 
systems (Fig. [2]^), we find that when increasing the cou- 
pling strength pxY beyond about 0.06, the fourth-largest 
Lyapunov exponent A4 gets negative. At about the same 
value, CPR starts increasing indicating the possible onset 
of PS. At fj,xY ~ 0.17, the third-largest Lyapunov expo- 
nent A3 becomes significantly negative, which is commonly 
interpreted as an indicator for the onset of GS. At the 
same time, we observe a sharp rise of JPR towards values 
of about 0.8, whereas Qj- immediately rises to values close 
to 1 (in fact, there is some weak overshooting because the 
edge densities of RNs and JRN differ necessarily, resulting 
in slightly different estimates of the corresponding network 
transitivities). This finding suggests that Qj- has better 
capabilities for detecting GS than JPR. 



The results are qualitatively similar when studying the 
same coupling configuration with both systems being in 
the non-coherent funnel regime (Fig. [2j3). Here, A4 be- 
comes negative already at rather low coupling strengths, 
whereas A3 < for fXxY ~ 0.3. The corresponding 
emergence of GS is accompanied by a continuous tran- 
sition of JPR and Qj- to values of about 1, where the 
JRN-based index again shows better convergence. More- 
over, for symmetric bidirectionally coupled funnel systems 
{fJ-XY — Pyx = P, Fig. [sjB), we observe qualitatively the 
same behaviour as in the unidirectional case, with the cor- 
responding transitions taking place at smaller values of the 
coupling strength (i.e., A4 < for /i > 0.02 and A3 < 
for ?i O-l''' ll''' )• Notably, the transition to GS as un- 
veiled by both JPR and Qj- is even more continuous than 
for unidirectional coupling, and the residual difference of 
JRN to the ideal value of 1 close to the transition point is 
even larger, whereas Qj- is at the same time already very 
close to 1. 

Studying two symmetrically coupled phase-coherent 
Rossler systems (Fig. |3]4_), the synchronisation scenario 
gets more complex due to the emergence of periodic win- 
dows as p is increased. Specifically, we find a first very 
small periodic window (Ai = 0, A2 < 0) close to the 
onset of PS at about p = 0.039, where CPR shows a 
sharp rise to values close to 1 17 . Note that in this 



as well as other supposedly periodic windows, we observe 
Tx ,Jy ,Tj ~ 0.75 consistent with the expected attrac- 
tor dimension of 1. At p > 0.04, also JPR and Qj- in- 
crease considerably, reaching plateau values in the phase- 
synchronised chaotic regime flTl. For p € [0.08,0.092], we 
find another more extended window of apparently periodic 
dynamics A2 < that probably corresponds to a regime 
of intermittent lag synchronisation 17 , whereas for even 



larger coupling strength, we infer the presence of GS since 
A3 < 0, JPR w 1 and Qr « 1. 

In order to illustrate the limitations of the proposed ap- 
proach as well as the established JPR index, we finally 
study the signatures of synchronisation for the mixed case 
of a phase-coherent Rossler system X coupled to a funnel 
oscillator Y, which can be taken as an example for syn- 
chronisation between two strongly non-identical systems. 
In the case of unidirectional coupling (Fig. ^p), we find 
A4 < for pxY ^0.1 accompanied by an increase of 
CPR indicating the presence of PS. At pxY ~ 0.3, also 
A3 approaches negative values, whereas both JPR and 
Qj- increase gradually, but do not approach values close 
to 1 even for relatively high coupling strengths. Thus, the 
three possible GS indicators A3 , JPR and Qj- do not give 
consistent results in this case, which implies that we can- 
not unambigously conclude the possible presence of GS 
here. A similar statement holds for the bidirectional case 
(Fig. |3p), which displays again a more complicated se- 
quence of transitions. 

Conclusions. — 25 years after the introduction of re- 
currence plots by Eckmann et al. 33 , the development 
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I^XY I^XT f^XY 



Figure 2: Results of synchronisation analysis for two unidirectionally {X — > Y) coupled Rossler systems being (A) both in the 
phase-coherent regime, (B) both in the funnel regime, and (C) in phase-coherent {X) and funnel regime (Y): the four largest 
Lyapunov exponents Ai, . . . , A4 estimated using the Wolf algorithm ^46j (using TV = 9, 000, 000 data points from simulations 
with step size h — 0.001 starting at T = 1,000); recurrence-based synchronisation indices CPR (black) and J PR (grey) |17| ; 
transitivities of individual and joint recurrence networks Tx (dark grey), Ty (light grey) and 7j (black); transitivity ratio Qj- 
(from top to bottom). 



of recurrence-based techniques still continues. With this 
paper, we would like to honour the seminal work by these 
authors and show how the paradigms of recurrence plots 
and complex networks have meanwhile been combined to 
pave the way for new concepts in complex systems analy- 
sis. 

We have introduced a new fully geometric approach 
for detecting GS between two coupled chaotic oscillators 
based on time series data. Specifically, we have used a 
complex network interpretation of the recurrence as well 
as joint recurrence matrices, the transitivity properties of 
which allow tracing changes in the effective dimensionality 
of the individual system's attractors as well as the attrac- 
tor of the combined higher-dimensional system. We have 
demonstrated that in the presence of GS, both systems 
effectively behave as one, with the dimensionality of the 
composed system approaching that of the individual ones. 

For the case of structurally similar systems such as 
slightly detuned chaotic oscillators with otherwise equal 
characteristic parameters, we have shown that the ratio 
Qj- between the JRN transitivity and the arithmetic mean 
of the individual systems' RNs transitivities quickly ap- 
proaches a value of 1 as GS takes place. In comparison 
with the conceptually related (i.e., recurrence plot-based) 
measure JPR, our geometric characteristic Qj- displays 
a better convergence towards the value of 1 expected for 
GS. In turn, for structurally different systems, both indices 
fail to clearly indicate the transition to GS as unveiled by 



studying the Lyapunov spectrum. 

In general, the proposed method has the advantage of 
being applicable to relatively short time series of length 
N < O(IO^), whereas the computational requirements 
for numerically estimating the Lyapunov spectrum are 
far higher and usually not met in case of real-world ap- 
plications. In this spirit, we conclude that the proposed 
method has great potentials for future applications to ob- 
servational data from different fields of research. 

^ ^ ^ 
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